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SECTION  I 
INTRODUCTION 

In  the  contemporary  technology  of  structural  composite  materials, 
major  deficiencies  exist  with  respect  to  our  ability  to  determine  the  stress 
field  within  a multilayered  composite  laminate.  In  most  cases,  even  a 
superficial  or  qualitative  understanding  of  the  nature  of  the  stress  field  in 
regions  of  steep  stress  gradients  has  not  been  established.  A notable  ex- 
ception is  the  approximate  treatment  [1,2]  of  the  classical  free  edge  prob- 
lem in  laminate  elasticity  [3,4],  however,  similar  treatments  have  not 
been  advanced  for  other  laminate  stress  concentration  problems.  In  the 
absence  of  a practical  means  of  laminate  stress  analysis,  it  is  not  possible 
to  develpp  an  understanding  and  general  characterization  of  the  various 
insidious  failure  modes  which  have  been  demonstrated  in  composite  laminates 
[2,  5,  6, 7, 8].  The  latter  references  all  document  heterogeneous  damage 
development  which  varies  through  the  laminate  thickness  and  suggest  the 
importance  of  defining  the  stress  field  within  each  layer,  particularly  in 
regions  of  stress  concentration,  where  classical  lamination  theory  [9,  10] 
errs  badly  [4,11,12].  This  is  contrasted  to  the  ultimate  failures  of  cer- 
tain laminates  in  the  presence  of  stress  risers,  which  only  seem  to  depend 
upon  the  overall  laminate  properties  [ 13,  14]  , rather  than  the  details  of 
the  stacking  sequence,  although  the  details  of  the  damage  development  and 
growth  in  these  laminates  are  a function  of  stresses  in  the  individual  layers 
[15]. 

Although  the  finite  element  method  is  widely  used  in  the  design  of 
practical  composite  structural  elements,  e.g.  [16],  its  application  is 
limited  to  determination  of  force  (per  unit  length),  rather  than  stress  dis- 
tributions, This  is  accomplished  through  the  assumption  of  a simplified 
displacement  field-that  which  is  assumed  in  classical  lamination  theory. 

This  permits  one  to  define  effective  elastic  properties  of  the  laminate  as  a 
whole,  and  to  compute  integrated  values  of  the  in-plane  stress  components 
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across  the  laminate  thickness.  In  fact,  once  the  force  distribution  has 
been  defined,  the  associated  layer  stresses  may  be  computed,  however,  as 
mentioned  earlier,  this  procedure  is  not  generally  reliable.  At  present, 
therefore,  determinations  of  failure  loading  and  mode  of  failure  in  practi- 
cal composite  structures  are  based  primarily  on  experimentation  with  pro- 
totype bodies.  These  comments  are  not  made  with  the  intention  of  dis- 
paraging the  contemporary  practice  since  the  presence  of  structural  com- 
plexity in  the  form  of  holes,  connections,  edges,  and  discontinuities  in 
thickness  of  many-layered  composites  may  overwhelm  any  attempted  analy- 
tical technique. 

The  basic  limitation  imposed  by  numerical  solutions  to  the  laminate 
elasticity  problem  has  been  illustrated  in  a recent  paper  by  Wang  and 
Crossman  [17]  in  their  treatment  of  the  free  edge  class  of  boundary  value 
problems.  In  order  to  achieve  a realistic  prediction  of  the  stress  field  by 
use  of  the  finite  element  method,  sixteen  elements  in  the  thickness  direc- 
tion were  required  within  each  layer  in  the  region  of  steep  stress  gradients. 
A total  of  196  elements  per  layer  were  employed.  To  accomodate  such  a 
large  array,  it  became  necessary  to  employ  a special  matrix  storage 
scheme  for  the  purpose  of  reducing  computer  storage  and  running  time. 
Similarly,  Pipes  [18]  required  a lengthy  extrapolation  procedure  in  con- 
junction with  the  finite  difference  method  in  order  to  achieve  a satisfactory 
solution  for  a two  layer  free  edge  boundary  value  problem. 

Another  approach  to  the  problem  of  laminate  stress  analysis,  e.g., 
Rybicki  [19],  utilizes  larger  elements  that  possess  a more  complex  stress 
field.  Stanton,  Crane,  and  Neu  [20]  employ  a tricubic  isoparametric 
discrete  element  and  a system  to  automate  the  construction  of  finite  ele- 
ment models.  The  latter  approach  effects  an  appreciable  reduction  in  data 
input  requirements.  This  added  efficiency  seems  to  be  the  major  advantage 
of  "large  element"  formulations,  as  the  number  of  degrees  of  freedom  is 
comparable  to  those  employed  in  the  more  routine  methods.  Despite  the 
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refinements,  however,  computed  laminate  stress  fields  are  not  universally 
accurate,  as  another  feature  of  laminate  elastic  analysis,  i.e.,  the  presence 
of  stress  singularities,  is  a severe  obstacle  to  the  execution  of  numerical 
elasticity  solutions.  Unfortunately,  the  order  of  these  singularities  has  not 
been  defined  for  anisotropic  layers,  i.e.,  the  extension  of  Bogy's  work  on 
isotropic  quarter-planes  [21]  has  not  been  accomplished. 

On  the  other  hand,  approximate  theories  have  been  proposed  in 
attempts  to  execute  realistic  laminate  stress  analysis.  The  most  popular 
of  these  is  the  aforementioned  classical  lamination  theory,  op.  cit.  [9.  10]  . 
While  this  theory  has  been  shown  to  yield  reasonably  accurate  stress  calcu- 
lations in  certain  classes  of  boundary  value  problems  under  a limited  range 
of  geometric  parameters  [22,23,24],  its  assumptions  are  too  restrictive 
for  general  application.  Theories  incorporating  laminate  "shear  deforma- 
tion" [25,26]  lead  to  accurate  deflections  in  bending  problems  but  offer  no 
hope  for  improved  stress  computations,  c.f.  [26].  The  higher-order  plate 
theory  derived  by  Whitney  and  Sun  [27]  was  applied  by  Pagano  [28]  to 
examine  the  interlaminar  normal  stress  distribution  in  the  free  edge  bound- 
ary value  problem  but  only  on  a plane  of  symmetry. 

The  prominent  common  feature  of  the  theories  discussed  in  the  pre- 
vious paragraph  is  an  assumed  displacement  field  that  is  continuous  across 
the  entire  laminate  thickness.  The  theories  differ  only  in  their  specific 
choice  of  the  assumed  displacements.  This  displacement  assumption,  how- 
ever, guarantees  discontinuous  tractions  at  interfaces  between  layers  of 
different  elastic  moduli  except  under  elementary  loading  conditions.  Fur- 
ther, the  edge  (traction)  boundary  conditions  appropriate  to  this  approach 
are,  in  general,  insufficient  to  guarantee  equilibrium  of  sub- regions  con- 
taining the  edge  boundary  under  the  known  (pointwise)  tractions,  op.  cit.  [28]=! 

=f:In  order  to  guarantee  equilibrium  of  a given  sub- region  containing  the  edge 
boundary,  we  must  have  the  freedom  to  prescribe  at  least  5 traction  bound- 
ary conditions  (3  force  components  and  2 couples)  on  its  edge.  Therefore, 
if  we  wish  to  guarantee  equilibrium  of  each  layer  of  a laminate,  we  need  at 
least  5N  edge  boundary  conditions,  where  N is  the  number  of  layers  in  the 
laminate. 
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Hence,  use  of  this  displacement  approach  and  possible  extensions  to  allow 
even  higher  order  variations  through  the  thickness,  is  unacceptable  for 
laminate  stress  field  analysis. 

Another  class  of  approximate  laminate  theories  represent  attempts  to 
generalize  those  discussed  above  and  are  based  upon  the  assumption  that  the 
displacement  components  are  linear  functions  of  z,  the  thickness  coordinate, 
within  each  layer.  In  this  formalism  then,  the  displacements  are  piecewise 
continuous  functions.  Among  the  theories  which  rely  on  this  approach  are 
the  so-called  effective  stiffness  theories  pioneered  by  Sun,  Achenbach,  and 
Herrmann  [29,30]-  Sun  and  Whitney  [31]  treated  various  theories  in  this 
class  and  demonstrated  that,  given  displacement  continuity  at  all  interfaces, 
the  number  of  field  equations  depend  on  N (number  of  layers)  only  when 
traction  continuity  at  interfaces  is  ignored.  Otherwise,  the  number  of  field 
equations  is  only  dependent  on  the  generality  of  the  initial  assumption,  i.e., 
whether  the  linear  term  in  z for  transverse  displacement  w is  included  or 
dropped.  Hence  the  number  of  field  equations  is  constant  for  add  laminates. 
Since  the  same  statement  can  be  made  with  respect  to  the  number  of  edge 
boundary  conditions,  the  deficiency  of  the  aforementioned  theories  [9,  10, 
25,26,27]  with  respect  to  sub-region  equilibrium  applies  to  the  present 
class  as  well.  The  latter  theories  do,  however,  yield  a more  realistic  de- 
termination of  effective  laminate  dispersion  characteristics,  which  provided 
the  motivation  for  their  development.  The  assumption  of  piecewise  linear 
displacements,  as  well  as  w = w(x,y),  leads  to  the  theory  of  Srinivas  [32], 
in  which  the  number  of  field  equations  and  edge  boundary  conditions  does 
depend  upon  the  number  of  layers.  Therefore,  interface  traction  continuity 
conditions  cannot  be  satisfied  [31].  Furthermore,  in  this  theory,  the  inter- 
laminar normal  stress,  which  has  been  shown  to  be  responsible  for  delam- 
ination failures  [2],  has  been  neglected.  Despite  the  accurate  results 
obtained  for  vibration  frequencies,  deflection,  and  axial  stress  in  the  parti- 
cular bending  problems  treated  in  [32],  the  theory  is  not  generally  applicable 
for  laminate  stress  analysis. 
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Finally,  an  approach  suggested  by  Pagano  [28],  i.e.,  treatment  of 
each  layer  as  a plate  governed  by  the  Whitney- Sun  theory  [27],  lacks 
generality  since  it  can  be  shown  that  less  than  5N  edge  traction  boundary 
conditions  are  available  in  that  approach.  Furthermore,  as  a consequence 
of  interface  displacement  continuity,  the  natural  edge  traction  boundary 
conditions  are  coupled,  i.e.,  they  involve  functions  of  the  tractions  acting 
on  two  (adjacent)  layers.  Thus,  extension  of  this  approach  by  allowing 
higher  order  displacements  is  not  advisable,  since,  as  in  the  previous 
approaches,  proper  equilibrium  of  each  layer  under  its  prescribed  tractions 
cannot  be  enforced. 

The  previous  discussion  has  defined  a clear  need  to  examine  new 
approaches  for  laminate  stress  analysis.  All  known  approximate  laminate 
theories  are  based  upon  assumed  displacement  fields,  which  as  we  have 
seen,  lead  to  results  lacking  credibility.  Therefore,  in  this  work,  we  shall 
set  down  requirements  to  be  satisfied  by  an  acceptable  laminate  field 
theory  and  proceed  to  develop  a self-consistent  theory  in  accord  with  the 
requirements,  which  are:  a)  All  six  stress  components  are  non- zero  in 
general,  b)  Traction  and  displacement  continuity  conditions  at  interfaces 
between  adjacent  layers  are  satisfied,  c)  Consider  a region  within  the  lamin- 
ate that  is  arbitrarily  located  except  that  it  is  bounded  by  any  two  of  the 
parallel  interfaces  of  the  laminate.  We  shall  require  that  the  computed 
stress  field  acting  on  the  surfaces  of  the  arbitrary  region,  in  conjunction 
with  the  prescribed  traction  boundary  conditions  (pointwise,  in  the  elasticity 
sense)  on  those  portions  of  the  external  laminate  boundary  which  lie  in  the 
region,  satisfy  the  conditions  of  vanishing  resultant  force  and  moment  iden- 
tically. Thus,  every  layer  must  satisfy  this  requirement,  which  we  shall 
henceforth  refer  to  as  "layer  equilibrium".  This  implies  that  the  appro- 
priate force  variables  in  the  field  theory  are  force  and  moment  resultants 
(per  unit  length)  acting  on  the  cross  sections  of  a layer  and  interlaminar 
stresses  on  its  interfacial  surfaces.  Although  sub- regions  not  bounded  by 
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interfacial  planes  need  not  satisfy  layer  equilibrium,  in  problem  solving, 
additional  interfaces  may  be  introduced  conceptually  to  improve  solution 
accuracy.  In  fact,  we  may  view  the  purpose  of  this  work  as  an  examination 
of  effectiveness  of  mathematical  laminate  models  in  which  the  response  is 
defined  in  terms  of  force  and  moment  resultants  and  interlaminar  stresses. 

Although  the  above  requirements  do  not  define  a unique  theory,  we 
shall  treat  the  simplest  theory  within  this  class  in  the  present  work.  The 
theory  is  based  upon  a variational  theorem  derived  by  Reissner  [33J  and 
permits  the  treatment  of  discontinuous  interfaces,  i.e.,  interface  cracks. 
Known  solutions  for  the  free  edge  boundary  value  problem  in  laminate 
elasticity,  where  pronounced  stress  gradients  occur,  will  be  utilized  to 
assess  the  consequences  of  the  present  concepts. 
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SECTION  II 

VARIATIONAL  PRINCIPLE  FOR  LAMINATES 

The  physical  problem  of  interest  in  the  present  work  is  that  of  a lam- 
inate which  is  built  of  anisotropic  elastic  layers  of  uniform  thickness  and  is 
subjected  to  prescribed  tractions  and/or  displacements  on  its  boundary. 

The  body  is  bounded  by  a cylindrical  edge  surface  and  upper  and  lower  faces 
that  are  parallel  to  the  interfacial  planes.  Since  it  is  necessary  to  consider 
both  traction  and  displacement  continuity  conditions  at  the  various  inter- 
faces, it  is  logical  to  examine  Reissner's  variational  theorem  [33]  as  a 
mechanism  to  develop  the  appropriate  field  equations. 

Reissner  has  shown  that  the  governing  equations  of  elasticity  can  be 
obtained  as  a consequence  of  the  variational  equation 

6J  = 0 (1) 

where 

J = / FdV  - / r.u.ds  (2) 


and 


F = ■r(r.,(u.  . + u.  .)  - W 

2 ij  i,j  3,1 


(3) 


In  the  above  equations,  W is  the  strain  energy  density  expressed  in  terms  of 
the  stresses  (r..(i,j  = 1,2,3),  V is  the  volume,  S the  entire  surface,  T.  the 
prescribed  tractions,  u^^  the  displacement  components,  and  S'  is  the  portion 
of  the  boundary  on  which  one  or  more  traction  components  are  prescribed. 

In  what  follows,  we  shall  letS"  represent  the  portion  of  the  boundary  on 
which  one  or  more  displacement  components  are  prescribed.  Summation 
over  the  range  of  repeated  subscripts  will  be  understood  in  this  work.  It  is 
also  understood  that  both  stresses  and  displacements  are  subjected  to  var- 
iation in  the  application  of  eq.  (1). 
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We  shall  now  express  the  form  of  eq.  (1)  for  the  laminated  body, 
which  is  composed  of  N layers,  the  volumes  of  which  are  represented  by 
Vj^(k  = 1,2,  • • • N).  For  conceptual  purposes,  we  may  let  the  layers  be 
numbered  consecutively  from  the  bottom  (k  = 1)  to  the  top  (k  = N).  Thus 
we  get,  by  definition 


J = 


IJ[^ 


IJ  1.3 


+ u. 


-|(k)  ^ 

W dV^  - j r.u.ds 


(4) 


where  the  superscript  (k)  attached  to  the  bracket  signifies  that  each  variable 
within  the  bracket  is  associated  with  the  kth  layer.  We  shall  also  incor- 
porate expansional  strains  e_  (34],  or  strains  produced  in  the  absence  of 
stress,  in  the  present  theory,  so  that 

W = W (a..,  e..)  (5) 

ij  13 


Substituting  (4)  into  (1),  making  some  trivial  manipulations,  and  applying 
the  Green- Gauss  Theorem  gives 


(k) 

where  S is  the  surface  enclosing  V,  and  r.  are  the  tractions  components 
k K 1 

acting  We  should  recognize  that  the  surfaces  and  contain 

a common  region,  namely,  the  interface  between  the  respective  layers. 
Hence,  we  shall  define  surface  as  the  portion  of  which  contains  the 
top  of  the  kth  layer.  We  also  let  represent  the  portion  of  1^  that  does 
not  belong  to  either  S'  orS".  Observing  that  the  edges  of  the  layers,  as 
well  as  the  top  of  the  Nth  layer  and  bottom  of  the  first  layer,  all  belong 
to  S,  eq.  ( 6)  may  be  expressed  as 


1. 


H 


[f 
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I /[(H^  -le) 

— 1 V,  -S' 


r N- 1 r 

I r.  6u.dS  + ^ / ( 

/.  ‘ ‘ ^ = ‘t 


(k)  (k)  (k+1)  (k+1) 

r.  6u.  + r.  6u.  ) 

.11  1 1 


Clearly,  the  vanishing  of  the  volume  integrals  requires  satisfaction  of  the 
equilibrium  equations  and  stress-displacement  relations  within  each  layer. 
The  vanishing  of  the  surface  integrals  on  S'  and  S"  require  that  one  term 
of  each  of  the  products  ( r u prescribed  at  each  point  on  S 

since  6Uj^  is  arbitrary  on  S'  and  it  vanishes  onS".  Finally,  the  integrals 
over  (k=l,2,"**N-  1)  vanish  when  tractions  and  displacements  are 
continuous  in  these  regions.  Hence,  eq.  (7),  which  represents  the  state- 
ment of  Reissner's  theorem  for  laminated  bodies,  will  be  applied  in  the 
derivation  of  our  approximate  laminate  theory  in  the  next  section. 
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SECTION  III 

DEVELOPMENT  OF  THEORY 


Consider  a single  layer  of  thickness  h within  the  laminate.  We  let  x 

and  y represent  the  coordinates  in  the  midplane  of  the  layer,  which  is 

bounded  by  the  planes  z = ± h/2  and  the  cylindrical  edge  surface  whose 

intersection  with  the  midplane  is  called  L.  The  region  enclosed  by  L will 

be  denoted  by  R.  The  interlaminar  stresses  u , r , and  r at  the  top  of 

z xz  yz  ^ 

the  layer  are  denoted  by  p^,  t^,  and  s^,  respectively,  while  the  corres- 
ponding stresses  at  the  bottom  of  the  layer  are  designated  as  p^^,  t^,  and  s^. 
Superscripts  (k),  which  identify  the  layers,  will  be  dropped  except  when 
they  are  needed  for  clarity. 


The  simplest  assumption  consistent  with  realistic  stress  analysis  is 

that  the  in-plane  stress  components*  are  linear  functions  of  z,  viz, 

N 


1 


12M  z 

X X 

h ,3 


N 12M  z 

+ 


(8) 


= r 


xy 


Nx 

h 


12M  z 



where  functions  of  x and  y only.  Obviously,  these  functions 

represent  the  usual  force  and  moment  resultants  arising  in  plate  theory. 

We  have  also  indicated  the  symbols  for  the  stress  components  in  contracted 
notation  since  this  system  will  be  convenient  for  future  developments. 


<1 


* Note  that  we  refrain  from  assuming  the  form  of  the  displacement  field 
in  accordance  with  the  objectionable  features  of  that  approach  described 
earlier. 
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We  now  substitute  (8),  along  with  the  values  of  the  interlaminar  stresses 
at  z = ± h/2,  into  the  differential  equations  of  equilibrium,  which  leads  to 
the  following  distributions 


(P1+P2) 

/I2z^  \ 

(Po-P,) 

(40z^ 

■— 

3N 

z 

N 

ro 

O'-*  - O'  = .4 

3 z 4 

4 I 

h/ 

2h 

V u2  / 

V h ' 

V h 

\ h / 

( s , + S-) 

/ 2 

3V 

^ 2 

z . 1 2' 

■ /l7  — 

A + 

y 

4z  \ 

4 yz  Z 1 

h 4 

V ■ 

7 ' 

2h 

( 

/2z  8z^ \ 

V^'h^  / 


i , z (S+V 

/12z^  1 

\ 3V 

o’_=r  . 

. 5 xz  ' 2 r h 4 

j 2h 

( ) 

where  the  shear  resultants  V , V and  the  functions  N , M 

X y z z 

given  by 


(N  , M ) 
z z 


(7  ( 1 , z)dz 
z 


are  functions  of  x and  y atone.  The  functions  on  the  right  hand  side  of 
eqs.  (8)  and  (9)  are  not  all  independent  as  they  are  related  via  equilibrium 
and  continuity  considerations,  but  these  relations  will  be  subsequently 
developed  by  means  of  the  variational  equation  (7). 

In  general,  the  strain  energy  density  of  an  elastic  anisotropic  body 
is  given  by 

W = 4s..(r.(r.  + (T.e.  (i,j  = 1,2’*‘6)  (11 

2 ij  1 j 11 

where  contracted  notation  has  been  employed,  with  S^.  representing  the 
compliance  matrix  and  e.  the  engineering  expansional  strain  components. 
Since  structural  composite  laminates  are  generally  built  such  that  each 
layer  possesses  a plane  of  elastic  symmetry  parallel  to  xy,  we  shall  treat 
this  material  class  (monoclinic)  in  this  work  although  generally  anisotropic 


m 


1 

M 


! 


i 


AFML-TR-77-n4 

layers  may  be  treated  without  difficulty.  For  monoclinic  symmetry  with 
respect  to  the  xy  plane,  the  compliance  matrix  takes  the  form 


'"ll 

"l2 

®22 

s..  = 

"l3 

"23 

"33 

SYMM. 

0 

0 

0 

"44 

(IZ’i 

0 

0 

0 

"45 

"55 

"26 

"36 

0 

0 

"66_ 

while 

^^23  = 

0 

(13) 

®5 

^^3  = 

0 

for  monoclinic  symmetry. 


We  now  substitute  eqs.  (8),  (9),  and  (11),  taking  account  of  ( 1 Z)  and 
(13),  into  the  variational  equation  (7),  which  yields 
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(lli  cont'd) 
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where  u,  v,  and  w are  the  x,  y,  and  z components  of  displacement,  re- 
spectively, e^,  e^,  e^,  and  e^^  are  the  engineering  expansional  strain 
components,  and  h represents  the  thickness  of  the  kth  layer. 

Since  we  have  made  no  assumption  regarding  the  displacement  field, 
the  integrations  with  respect  to  z in  (14)  give  rise  to  weighted  average  dis- 
placements and  displacements  at  the  interfaces.  Therefore,  we  make  the 
definitions 


(f,  f*.  f)  = 


f (1 

^ h^ 


where  f may  represent  either  u,  v,  or  w.  We  also  let  u^,  v^,  w^  represent 
the  displacement  components  at  the  top  of  the  layer  and  u^,  v^ , w^  the 
corresponding  functions  at  the  bottom  of  the  layer.  Furthermore,  for 
internal  consistency  in  the  theory,  we  express  the  prescribed  tractions  on 
the  appropriate  portions  of  boundary  L as  follows 

, / 12M  z\ 


, / 12M  z 

i fj 

h \ n,  ^2 


- '^2- 


where  n and  s are  local  coordinates,  which  are  respectively  normal  and 

tangent  to  L.  We  note  that  T and  T-  give  the  values  of  shear  stress  T at 

I ^ nz 

the  bottom  and  top  of  the  layer,  respectively.  No  restrictions  are  placed 
on  the  nature  of  the  boundary  tractions  and/or  displacements  over  the  re- 
mainder of  the  laminate  boundary.  We  may  now  perform  the  integration  of 
(14)  with  respect  to  z.  Taking  (15)  and  (16)  into  account,  we  get 


AFML-TR-77-n4 


i If 

k=i  “'y 


R 


6N  h _ 

"T  ("7  - he  - S N - S N - S N - S N )+  • ^ v, 

h 2 X X llx  12y  13z  16  xy  h 2 y 


6N 


-he  - S,_N  -S,,N  - S,,N  - S,,N  .)  + 


y “I2‘x  “22‘  y “23^  z ^26‘V'  h 


3w'”-  he  - S N - S N 
z 1 3 X Z 3 y 


"33" 


tS_,n  -s  n +— — (p,+p, 

5 33  z 36  xy  10  1 2 


6N 


(u,  + V,  )-  he  - s N 
2 y X xy  16  X 


- S- ,N  - S,,N  - S,  ,N 
26  y 36  z 66  xy 


126M^  ^2 

+ :r^  (“T 

,3  4 X llx  I2y  13z 


16  xy 
126M 


126M  ^2 

X , h_  =:= 


-SwM  )+  (— V-.  -S,,M  -S„M  - S,  M -S,,M  ) 


4 y IZ  X ZZ  y Z3  z Z6  xy 


I h (3<^-  w)-  S M - S,  M - ^ S,,M  - S,  ,M 
4 1 3 x 23  y 7 33  z 36  xy 


S,  h 

33 


(P,-P9) 


28  '^1  ^2 


12  6M 


XY 


— (u*.  +v*,  ) - S,,M  -S,,M 
4 y X 16  X 26  y 


(17 


So.M  - S,,M 
36  z 66  xy 


+ 6V 


3u=^  . 3 A 


+ 7 (w,  -w,  )-7r(S_V  +S^  V ) 
h 4 X X 5h  45  y 55  X 


S S 

""tI  + + V 


+ 6V 


3 '1'  3 / 

+7(w,  -w,  )--^(S.  V +S^  V ) 
h 4 y y 5h  44  y 45  X 


^44  ^45 


■|^(N  +N  +t--t,)--^(N  +N  +s^- 
2 x,x  xy,  y 2 1 2 xy,  x y,  y 2 


7 6u* 
h 


M +M  - V +-  (t,  +t,) 
x,x  xy,y  X 2 1 2 


r 


--6v 


M +M  - V +-^  c,  + (.  ) 
xy,x  y, y y 2 '1  2 


.'-■r 


'IW j* ' ' iijj 


AFML-TR-77-114 


(17  cont'd) 
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+ ( ^2-  "lU  Vl  + ^2-  —)  8 + E / t26u2+  6V2 
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(k)  (k)  (k+l)  (k+l)  (k+1)  (k+1)  (k+1)  (k+1) 

+ 6u^  - 6v^  - 5w^  dxdy  = 0 


(17  cont'd) 


where  the  region  refers  to  the  bottom  of  the  kth  layer,  while,  as  before, 
is  the  top  of  the  kth  layer,  and  again,  primes  refer  to  regions  in  which 
traction  boundary  conditions  are  prescribed  and  double  primes  denote  re- 
gions of  prescribed  displacement.  The  functions  which  are  integrated 
over  L'  and  L"  are  given  by  symbols  previously  defined,  except  that  in  this 
case,  n,  s,  and  z components  are  employed.  Also,  in  (17),  the  following 

contractions  are  introduced 
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The  appropriate  field  equations  and  boundary  conditions  in  the  present 

theory  follow  from  eq.  (17).  However,  further  manipulations  are  necessary 

to  develop  these  relations  as  all  variations  appearing  in  (17)  are  not  inde- 

(kT  (k)  (k) 

pendent.  In  particular,  the  interlaminar  stresses  t.,  s.,  and  p.  (i  = 1,2), 

(k)  (k)  (k)^  1 ^ 

as  well  as  the  interfacial  displacements  u^,  v^,  w^,  may  or  may  not  be 

independent  of  each  other,  depending  upon  the  nature  of  the  prescribed  con- 
ditions along  the  interfacial  planes.  To  examine  this  question,  consider  the 
contribution  to  the  variational  integral  5J  of  the  terms  involving  §p.  and 
5w.  (i  = 1,2)  at  the  interface  layers  k and  k+  I.  Let  this  portion  of  6J  be 
denoted  by  oj^  (k  = 1,2,  ••  • N-1).  The  following  treatment  can  be  applied 

in  the  same  manner  for  the  terms  involving  6t^  and  6u.  and  also  for  those 
involving  6Sj^  and  6v. . 


Reference  to  (17)  shows  that 
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we  shall  consider  in  detail  two  types  of  boundary  conditions  appropriate  to  a 
region  of  the  two  surfaces  constituting  an  interface,  viz.  , continuity  of  p and 
w,  and  prescribed  p on  one  surface  while  w is  specified  on  the  other.  From 
these  results,  the  appropriate  treatment  for  other  combinations  of  prescribed 
traction  and  displacement  will  be  quite  evident.  Considering  the  continuous 
interfacial  region  first,  we  have,  in  such  a region. 


19 
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(k+1) 

P, 


in  region  1"^^.  Substituting  eqs.  (ZO)  and  their  first  variations  into  (19) 
yields 

(k)  f r (k)  (k)  (k)  (k+1)  (k+1)  (k+1)  (k) 

K"  = JJ  ^^2  - ^33^2+  ^1-  ®33  (21) 


(k)_  (k) 

where  5J  represents  the  contribution  of  I to  the  integral  6J  . Since 
(k)  P k p 

Sp^  is  arbitrary  in  this  region,  the  function  in  parenthesis  must  vanish, 

hence  the  boundary  conditions  for  a continuous  interfacial  region  are  given 

by 

(k)  (k+l) 

P7  = P, 


(k)  (k)  (k) 

^2  ■ ^33  ^2 


(k+1)  (k+1)  (k+1) 
^1  - ^33  ^1 


We  observe  that  the  interfacial  displacements  w.  are  absent  from  the 
governing  equations  in  this  region. 

We  now  consider  an  interfacial  region  in  which  mixed  boundary  con- 
ditions are  prescribed,  e.g.. 


(k+1) 


(k+1) 
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(k+1)  (k)  (k)  (k+1) 

In  this  case,  6p  and  6w  are  arbitrary,  while  6p_  and  6w  vanish  in  the 

region.  Putting  these  results  into  (19)  while  dropping  integrals  which  do 

not  pertain  to  the  given  region  leads  to  the  expressions 


(k+1)  (k+1)  (k+1)  (k+1) 

- ®33  = ®1 


(k)  (k+1) 

which  characterize  a region  where  p^  and  w^  are  prescribed.  Equations 

analogous  to  (22)  and  (24)  may  be  derived  for  the  remaining  interlaminar 

stresses  t.,  s.  and  their  associated  displacements  u.,  v.,  respectively. 

Other  combinations  of  boundary  conditions  are  treated  in  similar  fashion. 

The  results  will  be  presented  later.  For  interface  conditions  of  the  type 

(23),  we  see  that  the  only  interfacial  displacement  that  will  appear  in  the 

(k+1) 

governing  equations  is  w^,  which  is  a prescribed  function.  It  follows  that 
this  is  true  in  general,  i.e.,  only  prescribed  interfacial  displacement 
functions  enter  the  governing  equations,  hence,  interfacial  displacements 
are  not  treated  as  dependent  variables  in  the  present  theory. 

Returning  to  eq.  (17),  we  now  develop  the  appropriate  field  equations 
and  boundary  conditions  by  setting  to  zero  the  coefficients  of  the  arbitrary 
functions  (first  variations).  The  field  equations,  which  consist  of  the  elastic 
constitutive  relations  and  the  differential  equations  of  equilibrium,  must  be 
satisfied  within  each  layer  and  are  given  by 
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(25  cont'd) 
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V + V + 
x.x  y,y 


+ 5(Pi-P2)-  y+s^  y)  = 0 

(26  cont'd) 


The  vanishing  of  the  remaining  arbitrary  functions  within  the  region,  after 
taking  possible  interdependence  into  account  as  noted  earlier,  leads  to 
expressions  that  must  be  satisfied  within  each  interfacial  plane.  These 
interface  equations  depend  on  the  nature  of  the  prescribed  conditions  on  the 
interfacial  planes  and  are  given  by 


Interface  Conditions: 
a)  Continuity  (k  = 1,2,  •••  N-1) 


(k) 

(k+l) 

h 

= 

S 

(k) 

(k+l) 

"2 

"1 

(k) 

(k+l) 

P2 

- 

Pi 

(27) 

(k) 

(k)  (k) 

(k) 

(k) 

(k+l) 

(k+l)  (k+l) 

(k+l)  (k+l) 

P4  - 

44  4 

- ^45 

“4 

- "44  Q4  - 

Qc  = 0 

45  5 

(k) 

(k)  (k) 

(k) 

(k) 

(k+l) 

(k+l)  (k+l) 

(k+l)  (k+l) 

P5  - 

45  4 

- "55 

T5  4 

“5 

- S..  Q.  - 

45  4 

S55  Q3  = 0 

(k) 

(k) 

(k) 

(k+l)  (k+l)  (k+l) 

^2 

■ ®33 

^2 

+ Vj  - S33  Rj  = 0 

Prescribed  Tractions  and/or 

Displacements  (k  = 1,2 

. . .N-1) 

(k) 

(k) 

(k) 

(k)  (k)  (k) 

(k)  (k) 

‘2  = 

^2 

or 

P5 

- ®45  ^4  - ®55 

^5  = -"2 

(k) 

(k) 

(k) 

(k)  (k)  (k) 

(k)  (k)  (28) 

"2  = 

^2 

or 

P4 

- S..  T . - S.. 
44  4 45 

T = - ^ 

5 2 
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(k)  (k) 

Pz  ^ ^2 

(k+1)  (k+1) 

(k+1)  (k+1) 
"l  = ^1 

(k+1)  (k+1) 

Pi  = Pi 


(k)  (k)  (k)  (k) 

^2  ■ ®33  ^2  ■ ^"2 


(k+1)  (k+1)  (k+1)  (k+1)  (k+1)  (k  + 1) 
“5  ■ ®45  ^4  ■ ^55  ^5  ""  “l 


(k+1)  (k+1)  (k+1)  (k+1)  (k+1)  (k  + 1) 

Q.  - S.  . Q.  - S Q = V 
4 44  4 45  5 1 


(k+l)(k+l)(k+l)  (k+1) 

- ®33  ‘'l  = *1 


(28  cont'd) 


where  eqs.  (28)  are  to  be  understood  in  the  sense  that,  at  each  interface, 
any  combination  which  contains  one  equation  from  each  line  can  be  used  to 
represent  the  interface  boundary  conditions  in  any  region  of  the  interface 
for  which  (27)  are  not  prescribed. 

Finally,  vanishing  of  the  exterior  surface  integrals  in  (17)  lead  to 
the  definition  of  the  natural  boundary  conditions  in  this  theory,  i.e.. 


Boundary  Conditions: 

a)  Edge  Surface 

For  the  edge  surface,  one  term  from  each  of  the  following  products 
must  be  prescribed  for  each  layer  (superscripts  k are  omitted) 

. /'^n  ^+^2\  _ 


N u , 

N 

u , M u"'' , M u"'' 

n n 

ns 

s n n ns  s 

/ 2V  \ 

1 n\  A 

>'2-  '1' 

W'% 

\Y^z-  —) 

b)  Top  Surface 

The  boundary  conditions  on  the  top  surface  are  the  same  as  the  first 
three  lines  of  (28)  with  k = N. 


c)  Bottom  Surface 

The  boundary  conditions  on  the  bottom  surface  are  the  same  as  the 
last  three  lines  of  (28)  with  k = 0. 
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This  completes  the  development  of  the  present  theory.  We  observe 
that  the  governing  equations,  (25)  - (28),  plus  the  boundary  conditions  on 
the  top  and  bottom  surfaces,  constitute  a system  of  23N  equations  in  terms 
of  a like  number  of  unknowns.  The  system  can  be  reduced  to  1 3N  equations 
by  solving  (25)  for  the  force  and  moment  resultants  in  terms  of  the  weighted 
displacement  functions  and  interlaminar  stress  components  and  substituting 
into  the  remaining  equations.  From  (29),  we  see  that  7N  edge  conditions 
are  required  in  this  theory.  In  the  event  that  only  edge  tractions  are  pre- 
scribed in  a given  boundary  value  problem,  these  7N  edge  functions  may 
be  taken  to  be  the  3N  force  resultants,  2N  moment  resultants  and  2N  inter- 
laminar shear  stresses  at  the  top  and  bottom  of  every  layer. 

Clearly,  the  requirements  established  in  the  introductory  section  are 
all  satisfied  by  the  present  theory,  in  particular,  the  principle  of  "layer 
equilibrium".  Furthermore,  the  generality  of  the  interface  conditions, 
(27)and  (28),  allow  for  the  presence  of  interfacial  cracks  in  the  treatment 
of  specific  boundary  value  problems.  Finally,  the  usual  (physically  mean- 
ingful) equations  of  equilibrium  are  represented  by  the  first,  second, 
fourth,  and  fifth  of  (26),  along  with  a linear  combination  of  the  third  and 
seventh  of (26). 
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SECTION  IV 

COMPARISON  WITH  FINITE  ELEMENT  RESULTS 

In  this  section  we  shall  relate  the  response  predicted  by  the  present 
theory  to  that  given  by  numerical  elasticity  solutions  for  several  problems 
of  practical  and  theoretical  interest.  We  shall  treat  the  class  of  boundary 
value  problems  known  as  the  free  edge  problem  in  which  a laminate  of  finite 
width  is  subjected  to  a uniform  axial  strain  = e[4]  . The  origin  of 
coordinates  is  located  at  the  center  of  the  laminate  and  each  layer  is  rein- 
forced by  a system  of  parallel  fibers  oriented  at  an  angle  6 with  the  x-axis 
as  shown  in  Fig.  1.  The  fibers  in  the  various  layers  all  lie  in  planes  para- 
llel to  xy,  and  the  laminate  is  symmetric,  i.e.,  0(z)  = 0(-z).  In  the  analy- 

sis of  the  stress  field,  which  is  only  a function  of  y and  z,  each  layer  is 
treated  as  a homogeneous,  anisotropic  body  represented  by  its  effective 
moduli  and  stresses  will  be  denoted  by  functions  of  the  form  f(y,  z) . 

Comprehensive  results  based  upon  the  finite  element  method  have 
recently  been  presented  by  Wang  and  Crossman  [17]  for  this  class  of 
boundary  value  problems  in  laminate  elasticity.  Hence,  that  work  will  be 
employed  here  to  compare  specific  results  given  by  the  present  theory. 

Two  particular  laminates;  [0,  90],  in  which  the  values  of  0 in  consecutive 
layers  are  0°,  90°,  90°,  0°,  and  [±  45],  in  which  the  orientations  are 
45°,  -45°,  -45°,  45°,  will  be  examined  in  this  study.  The  layers  are  of 
equal  thickness  h,  the  laminate  width  is  2b  = 16h,  and  the  moduli  in  the 
planes  of  elastic  symmetry  of  each  layer  are  given  by 


Ejj  = 20  X lO^psi  , E22  = = 2.1  X lO^psi 

^12  = °13  = ^23  = .85  x 10^1 
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where  1,  2,  and  3 refer  to  the  fiber,  transverse,  and  thickness  directions, 
respectively,  and  for  example,  is  the  Poisson  ratio  measuring  strain 

in  the  transverse  direction  due  to  uniaxial  tension  in  the  fiber  direction. 

In  Figs.  2-5,  we  compare  various  features  of  the  response  for  the 
[±45]  laminate  as  given  by  the  present  theory  [ 35]  and  the  finite  element 
solution  of  [ 17]  . The  values  of  N in  these  figures  correspond  to  the  number 
of  sub- layers  used  in  the  present  theory  to  model  one-half  of  the  laminate. 
Thus,  N = 6 implies  that  each  physical  layer  of  thickness  h in  the  body  has 
been  modeled  by  three  sub- layers,  each  of  thickness  h/3,  while  N = 2 indi- 
cates that  each  physical  layer  is  treated  as  a unit. 

In  Figs.  2 and  3 are  shown  the  distribution  of  a and  r , respectively, 

X xy 

along  the  width  of  the  laminate  at  the  center  of  the  top  (physical)  layer.  The 
functions  given  by  the  present  theory  were  computed  via  eqs.  (8).  The 
results  for  N = 6 and  the  finite  element  solution  are  nearly  coincident  for 
all  values  of  y,  while  the  N = 2 results  differ  by  only  a few  percent  in  the 
boundary  layer  region. 

Even  the  N = 2 result  agrees  quite  well  with  that  of  the  finite  element 

solution  for  the  width  distribution  of  r at  the  ±45°  interface.  However, 

xz 

a singularity  is  expected  at  this  level  at  the  free  edge  [4,17].  The  presence 
of  a singularity  introduces  some  ambiguity  in  the  finite  element  solution,  so 
that  the  curve  given  here  involves  some  judgement  in  the  interpretation  of 
the  numerical  results.  On  the  other  hand,  the  singularity  is  manifested  in 
a different  way  in  the  present  theory,  i.e. , the  stress  component  r at  the 
singular  point  grows  with  increasing  values  of  N.  Whether  a finite  limit  is 
approached  for  large  N has  not  been  established.  This  situation  is  similar 
to  the  rise  in  average  stress  in  the  element  adjacent  to  the  singular  point 
as  the  element  size  decreases  in  the  finite  element  method.  However,  the 
present  theory  contains  no  singularity  (see  [ 35]),  consequently  the  compu- 
ted stress  distribution  is  an  exact  solution  in  this  formulation.  The  growth 
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of  stress  component  r^^(b,  h)  with  N is  shown  in  Table  1.  The  result  for 
N = 3 represents  the  average  given  by  the  case  in  which  the  lower  layer  is 
represented  by  2 sub- layers  and  the  upper  layer  by  one,  and  the  opposite 
situation,  although  the  two  results  are  nearly  identical.  The  same  inter- 
pretation is  invoked  for  N = 5.  Unfortunately,  because  of  the  magnitudes  of 
the  numbers  involved  in  the  solution  approach  employed  in  [ 35] , values  of 
N larger  than  6 could  not  be  considered.  Clearly,  the  manner  in  which 
singular  behavior  is  portrayed  in  the  present  theory  needs  further  study. 

In  particular,  the  approach  by  which  one  correlates  the  analytical  results 
with  delamination  failure  tests  needs  consideration. 

Although  displacement  components  are  not  dependent  variables  in  the 
present  theory,  the  weighted  displacement  functions  can  be  used  to  approxi- 
mate them  with  the  aid  of  an  assumed  variation  within  each  layer.  For 
example,  if  we  assume  that  axial  displacement  u is  a linear  function  of  z 
within  each  layer,  by  use  of  (15)  we  can  show  that 


where  z^^  is  measured  from  a local  coordinate  system  at  the  center  of  the 
kth  layer.  Agreement  between  this  approach  and  the  finite  element  result 
for  axial  displacement  distribution  across  the  width  of  the  top  surface  is 
quite  good  as  shown  in  Fig.  5. 

In  Fig.  6,  the  distribution  of  along  the  width  direction  on  the  cen- 
tral plane  (z  = 0)  of  the  [0,90]  laminate  is  shown.  Clearly,  the  present 
theory  with  N = 6 agrees  quite  well  with  the  finite  element  result,  while  the 
N = 2 result  appears  accurate  except  in  a region  very  close  to  the  free  edge. 
Not  shown  on  the  figure  is  the  result  for  N = 4,  which  has  a very  slight  hump 
near  the  free  edge  and  attains  a maximum  value  close  to  the  N = 6 result. 

Figure  7 illustrates  the  width- dependence  of  at  the  interface  be- 
tween the  0°  and  90°  layers,  where  in  contrast  to  Fig.  6,  a singularity  is 
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expected  at  the  free  edge  owing  to  the  discontinuity  in  elastic  properties. 

The  finite  element  solution  gives  strong  evidence  of  the  singularity  since 
extreme  variability  occurs  in  the  neighborhood  of  (b,h).  Because  of  this, 
the  finite  element  results  are  somewhat  subjective  in  this  region.  Again, 
the  N = 6 result  is  closer  to  the  finite  element  curve  than  that  of  N = 2.  As 
before  (Fig.  4),  the  present  theory  yields  finite  maximum  stresses  which 
appear  to  grow  monotonically  with  increasing  N at  the  singular  point. 

Comparative  results  for  the  distribution  of  r at  the  0°  - 90°  inter- 

yz 

face  are  shown  in  Fig.  8.  The  present  theory  satisfies  the  traction- free 
boundary  condition,  however,  whether  the  finite  element  solution,  or  indeed, 
an  exact  elasticity  solution,  satisfies  this  condition  (see  [36])  is  not  known. 
However,  generally  reasonable  agreement  can  be  observed.  According  to 
the  new  theory,  it  appears  that  the  function  is  approaching  a finite  peak  value, 
although  we  cannot  be  certain  until  the  solution  for  larger  values  of  N is 
determined. 

Variation  of  transverse  displacement  v at  the  top  surface  is  shown  in 
Fig.  9.  The  values  in  the  present  theory  were  defined  by  approximating 
the  layer  displacement  as  a linear  function  of  z,  which  leads  to  an  equation 
of  the  same  form  as  (30).  Excellent  agreement  is  seen  to  occur  between 
the  two  solution  techniques. 

An  extremely  steep  stress  gradient  at  an  (apparent)  singularity  in 
(T^  was  reported  by  Rybicki  and  Pagano  [ 37]  for  a free  edge  problem  in 
which  one  layer  was  isotropic.  Using  the  moduli  given  in  [37],  compara- 
tive results  are  shown  for  the  distributions  of  <t ^ at  the  interface  (Fig.  10) 
and  central  plane  (Fig.  11),  where  comparable  agreement  with  previous 
results  can  be  observed. 


AFML-TR-77-n4 


I 


CONCLUDING  REMARKS 


We  have  derived  an  approximate  theory  for  the  stress  analysis  of 
laminated  bodies  which  resolves  the  difficulties  involved  in  previous  theories 
based  upon  assumed  displacement  fields.  This  theory  is  based  upon  Reissner's 
variational  principle  and  assumed  in-plane  stresses  that  are  linear  functions 
of  thickness  coordinate  z within  each  layer.  While  the  appearance  of  1 3N 
field  equations  and  7N  edge  conditions  may  seem  to  be  overly  cumbersome 
in  actual  problem  solving,  this  level  of  detail  is  required  to  compute  real- 
istic global  stress  fields.  The  present  theory  guarant<^es  satisfaction  of 
"layer  equilibrium"  and  allows  the  prescription  of  combinations  of  inter- 
facial  tractions  and  displacements  which  permit  treatment  of  such  condi- 
tions as  interfacial  continuity  or  cracks. 

Comparison  with  existing  solutions  of  the  laminate  free-edge  class  of 
boundary  value  problems,  in  which  very  steep  stress  gradients  occur,  has 
led  to  encouraging  results.  Although  certain  highly  localized  details  of  the 
stress  field  have  been  expunged  when  each  layer  was  modeled  as  a single 
unit,  this  approach  may  be  adequate  for  purposes  of  structural  design.  If 
this  is  not  the  case,  based  on  the  present  study,  the  introduction  of  two  or 
three  sub- layers  will  produce  dramatic  improvements  in  accuracy.  Alter- 
natively, one  may  incorporate  higher  order  terms  in  z into  eqs.  (8)  to 
develop  more  accurate  theories  satisfying  the  basic  requirements  set  forth 
here.  Such  theories  may  eliminate  the  need  for  the  use  of  sub-layers,  but 
will  obviously  lead  to  greater  complexity  in  the  solutions  of  specific  bound- 
ary value  problems. 

The  situation  regarding  singularities  remains  somewhat  nebulous 
since  the  precise  nature  of  the  singular  stress  field  in  the  vicinity  of  an 
interface  at  an  edge  has  not  been  established.  We  can  state  however,  that 
the  finite  element  solution  can  be  severely  hampered  by  the  presence  of 
elastic  stress  singularities,  and  stress  field  determination  in  their  vicinity 


30 


tains  no  edge  singularities,  an  advantage  from  the  problem  solving  view- 
point, however,  examples  have  demonstrated  a tendency  for  the  computed 
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maximum  stress  to  grow  with  decreasing  sub- layer  thickness.  Thus, 
problem  solving  has  become  simplified,  but  a method  to  interpret  the  stress 
predictions  needs  to  be  developed.  We  should  notice,  however,  that  the 
singularities  given  in  effective  modulus  theories  are  mathematical  artifacts 
in  the  treatment  of  fiber  reinforced  laminated  bodies.  This  has  been  dis- 
cussed in  [ 37,  38J , where  support  was  given  to  the  use  of  integrated  stresses 
rather  than  point  stresses  in  failure  analysis.  This  point,  together  with 
the  automatic  satisfaction  of  layer  equilibrium,  as  well  as  the  capability 
for  objective  determination  of  laminate  stress  fields,  favor  the  use  of  the 
present  theory  over  approaches  based  upon  numerical  solutions  of  the 
elasticity  equations.  Unfortunately,  this  work,  along  with  that  of  [ 35j  and 
[17],  demonstrate  the  extreme  difficulties  associated  with  attempts  to 
realistically  define  the  stress  fields  in  laminates  consisting  of  very  many 
layers . 
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Figure  1.  Laminate  Geometry 
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Distribution  of  r Along  Center  of  Top  Layer  (z  = 3^  h) 


Figure  4.  Distribution  of  t Along  45°  Interface  (z  - h) 
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Figure  11.  Distribution  of  <r.  Along  Central  Plane  (z  = 0) 


